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ABSTRACT 

In these lectures we introduce the linear operators of fractional integration and frac- 
tional differentiation in the framework of the Riemann-Liouville fractional calculus. 
Particular attention is devoted to the technique of Laplace transforms for treating 
these operators in a way accessible to applied scientists, avoiding unproductive gen- 
eralities and excessive mathematical rigor. By applying this technique we shall derive 
the analytical solutions of the most simple linear integral and differential equations of 
fractional order. We shall show the fundamental role of the Mittag-Lefffer function, 
whose properties are reported in an ad hoc Appendix. The topics discussed here 
will be: (a) essentials of Riemann-Liouville fractional calculus with basic formulas 
of Laplace transforms, (b) Abel type integral equations of first and second kind, (c) 
relaxation and oscillation type differential equations of fractional order. 
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1. INTRODUCTION TO FRACTIONAL CALCULUS 

1.1 Historical Foreword 

Fractional calculus is the field of mathematical analysis which deals with the 
investigation and applications of integrals and derivatives of arbitrary order. The 
term fractional is a misnomer, but it is retained following the prevailing use. 

The fractional calculus may be considered an old and yet novel topic. It is an 
old topic since, starting from some speculations of G.W. Leibniz (1695, 1697) and 
L. Euler (1730), it has been developed up to nowadays. A list of mathematicians, 
who have provided important contributions up to the middle of our century, includes 
P.S. Laplace (1812), J.B.J. Fourier (1822), N.H. Abel (1823-1826), J. Liouville (1832- 
1873), B. Riemann (1847), H. Holmgren (1865-67), A.K. Griinwald (1867-1872), A.V. 
Letnikov (1868-1872), H. Laurent (1884), RA. Nekrassov (1888), A. Krug (1890), J. 
Hadamard (1892), O. Heaviside (1892-1912), S. Pincherle (1902), G.H. Hardy and 
J.E. Littlewood (1917-1928), H. Weyl (1917), P. Levy (1923), A. Marchaud (1927), 
H.T. Davis (1924-1936), A. Zygmund (1935-1945), E.R. Love (1938-1996), A. Erdelyi 
(1939-1965), H. Kober (1940), D.V. Widder (1941), M. Riesz (1949). 

However, it may be considered a novel topic as well, since only from a little more 
than twenty years it has been object of specialized conferences and treatises. For the 
first conference the merit is ascribed to B. Ross who organized the First Conference 
on Fractional Calculus and its Applications at the University of New Haven in June 
1974, and edited the proceedings, see [1]. For the first monograph the merit is 
ascribed to K.B. Oldham and J. Spanier, see [2], who, after a joint collaboration 
started in 1968, published a book devoted to fractional calculus in 1974. Nowadays, 
the list of texts and proceedings devoted solely or partly to fractional calculus and its 
applications includes about a dozen of titles [1-14], among which the encyclopaedic 
treatise by Samko, Kilbas & Marichev [5] is the most prominent. Furthermore, we 
recall the attention to the treatises by Davis [15], Erdelyi [16], Gel'fand & Shilov 
[17], Djrbashian [18, 22], Caputo [19], Babenko [20], Gorenfio & Vessella [21], which 
contain a detailed analysis of some mathematical aspects and/or physical applications 
of fractional calculus, although without explicit mention in their titles. 

In recent years considerable interest in fractional calculus has been stimulated 
by the applications that this calculus finds in numerical analysis and different areas 
of physics and engineering, possibly including fractal phenomena. In this respect 
A. Carpinteri and F. Mainardi have edited the present book of lecture notes and 
entitled it as Fractals and Fractional Calculus in Continuum Mechanics. For the 
topic of fractional calculus, in addition to this joint article of introduction, we have 
contributed also with two single articles, one by Gorenfio [23], devoted to numerical 
methods, and one by Mainardi [24], concerning applications in mechanics. 
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1.2 The Fractional Integral 

According to the Riemann-Liouville approach to fractional calculus the notion of 
fractional integral of order a (a > 0) is a natural consequence of the well known 
formula (usually attributed to Cauchy), that reduces the calculation of the n— fold 
primitive of a function f(t) to a single integral of convolution type. In our notation 
the Cauchy formula reads 



where IN is the set of positive integers. From this definition we note that fn{t) 
vanishes at t = with its derivatives of order 1,2, . . . ,n ~ 1. For convention we 
require that f{t) and henceforth fn{t) be a causal function, i.e. identically vanishing 



In a natural way one is led to extend the above formula from positive integer 
values of the index to any positive real values by using the Gamma function. Indeed, 
noting that (n — 1)! = T{n) , and introducing the arbitrary positive real number a , 
one defines the Fractional Integral of order a > : 



where H is the set of positive real numbers. For complementation we define J := / 
(Identity operator), i.e. we mean J° f{t) = f{t) . Furthermore, by J"/(0"'') we mean 
the limit (if it exists) of J'^f{t) for t — > 0"'" ; this limit may be infinite. 
Remark 1 : 

Here, and in all our following treatment, the integrals are intended in the generalized 
Riemann sense, so that any function is required to be locally absolutely integrable 
in IR"*". However, we will not bother to give descriptions of sets of admissible func- 
tions and will not hesitate, when necessary, to use formal expressions with generalized 
functions (distributions), which, as far as possible, will be re- interpreted in the frame- 
work of classical functions. The reader interested in the strict mathematical rigor is 
referred to [5] , where the fractional calculus is treated in the framework of Lebesgue 
spaces of summable functions and Sobolev spaces of generalized functions. 
Remark 2 : 

In order to remain in accordance with the standard notation / for the Identity oper- 
ator we use the character J for the integral operator and its power J". If one likes 
to denote by I'^ the integral operators, he would adopt a different notation for the 
Identity, e.g. I, to avoid a possible confusion. 




(1.1) 



for t < . 




(1.2) 
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We note the semigroup property 

rjp = r+p^ tt,/3>o, (1.3) 

which imphes the commutative property — J'^J^, and the effect of our opera- 

tors J" on the power functions 

J°'t^= r!^^l^^ «>0, 7>-l, *>0. (1.4) 

r(7 + l + a) 

The properties (1.3-4) are of course a natural generahzation of those known when 
the order is a positive integer. The proofs, see e.g. [2], [5] or [10], are based on the 
properties of the two Eulerian integrals, i.e. the Gamma and Beta functions, 

/>oo 

T{z):^ e-'^u'-^du, T{z + 1) = zT{z) , Re{z}>0, (1.5) 
Jo 

B{p,q):= j\l-uY-'u'^-Uu = ^Ml^ ^ ^(g,^) , Re{p,g}>0. (1.6) 
It may be convenient to introduce the following causal function 

^a— 1 

^a{t):=^, «>0, (1.7) 

where the suffix + is just denoting that the function is vanishing for t < . Being 
q: > , this function turns out to be locally absolutely intcgrable in H"*" . Let us now 
recall the notion of Laplace convolution., i.e. the convolution integral with two causal 
functions, which reads in a standard notation f{t) * g{t) := f{t — r) g{T) dr = 
9{t) * fit) . 

Then we note from (1.2) and (1.7) that the fractional integral of order o: > can 
be considered as the Laplace convolution between $a(^) and f{t) , i.e. 

r f{t)=^^{t) * f{t), a>0. (1.8) 

Furthermore, based on the Eulerian integrals, one proves the composition rule 

$«(t) * ^/3(t) = (t), «,/3>0, (1.9) 

which can be used to re-obtain (1.3) and (1.4). 

Introducing the Laplace transform by the notation C {f{t)} := J^g~^* f{t) dt — 
f{s) , s e C , and using the sign to denote a Laplace transform pair, i.e. f{t)-^f{s) , 
we note the following rule for the Laplace transform of the fractional integral, 

^"/W-^, «>0, (1.10) 

which is the straightforward generalization of the case with an n-fold repeated integral 
[a = n). For the proof it is sufficient to recall the convolution theorem for Laplace 
transforms and note the pair ^a{t) -j- 1/s" , with a > , see e.g. Doetsch [25]. 
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1.3 The Fractional Derivative 

After the notion of fractional integral, that of fractional derivative of order a 
(a > 0) becomes a natural requirement and one is attempted to substitute a with 
—a in the above formulas. However, this generalization needs some care in order to 
guarantee the convergence of the integrals and preserve the well known properties of 
the ordinary derivative of integer order. 

Denoting by with n e IN , the operator of the derivative of order n , we first 
note that 

DT = I, rD^^^I, nGlN, (1.11) 

i.e. is left-inverse (and not right-inverse) to the corresponding integral operator 
. In fact we easily recognize from (1.1) that 



fc=0 



k\ 



(1.12) 



As a consequence we expect that D°' is defined as left- inverse to J°'. For this 
purpose, introducing the positive integer m such that m — 1 < a < m , one defines 
the Fractional Derivative of order a > : D'^ f{t) := jm-a j^^-j ^ namely 



D"f{t) -.^l 



df^ 
d"" 



1 r /M 

-a) Jo (t- t)«+i-^ 



r(m 

fit) 



m — 1 < a < m . 



(1.13) 



a = m . 



(1.14) 



(1.15) 



Defining for complementation = = I , then we easily recognize that 

D'^J'^ = I, a > , 

and 

rf7-h 1) 

r(7 + l-a) 

Of course, the properties (1.14-15) are a natural generalization of those known when 
the order is a positive integer. Since in (1.15) the argument of the Gamma function 
in the denominator can be negative, we need to consider the analytical continuation 
oiT{z) in (1.5) to the left half-plane, see e.g. Henrici [26]. 

Note the remarkable fact that the fractional derivative f is not zero for the 
constant function f{t) = 1 if o; ^ IN . In fact, (1.15) with 7 = teaches us that 

L»«l = -f^, a>0, t>0. (1.16) 
1 (1 — a) 

This, of course, is = for a e INT, due to the poles of the gamma function in the 
points 0, —1, —2, . . .. 
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We now observe that an alternative definition of fractional derivative, orig- 
inally introduced by Caputo [19], [27] in the late sixties and adopted by Ca- 
puto and Mainardi [28] in the framework of the theory of Linear Viscoelasticity 
(see a review in [24]), is the so-called Caputo Fractional Derivative of order a > : 

f{t) := J""-" D"" f{t) with m-l<a<m, namely 

1 /("^)(r) 



Tim -a) Jo t-T)«+i-"^ ' ' , , 

— f(t) , a = m . 

df^ ^ ' 



This definition is of course more restrictive than (1.13), in that requires the absolute 
integrability of the derivative of order m. Whenever we use the operator we 
(tacitly) assume that this condition is met. We easily recognize that in general 

f{t) := J"^-" f{t) J^-'^ /(t) := /(t) , (1.18) 

unless the function jit) along with its first m — 1 derivatives vanishes at t = 0+. In 
fact, assuming that the passage of the m-derivative under the integral is legitimate, 
one recognizes that, for m — 1 < a < m and t > , 

m — 1 j.k—a 

D'^ m = m + j2 r(k-a+i) ^^'^^^^^ ' ^^-^^^ 

and therefore, recalling the fractional derivative of the power functions (1.15), 

/(^) - E ^ /^'^ (0^) 1 = ^* /(^) • (1-20) 

The alternative definition (1.17) for the fractional derivative thus incorporates the 
initial values of the function and of its integer derivatives of lower order. The sub- 
traction of the Taylor polynomial of degree m — 1 at t = O"*" from f{t) means a sort of 
regularization of the fractional derivative. In particular, according to this definition, 
the relevant property for which the fractional derivative of a constant is still zero, i.e. 

D^1 = 0, a>0. (1.21) 

can be easily recognized. 

We now explore the most relevant difi^erences between the two fractional deriva- 
tives (1.13) and (1.17). We agree to denote (1.17) as the Caputo fractional derivative 
to distinguish it from the standard Riemann-Liouville fractional derivative (1.13). 
We observe, again by looking at (1.15), that 

j^a^a-l^Q^ q;>0, t>0. (1.22) 
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From (1.22) and (1.21) we thus recognize the following statements about functions 
which for t > admit the same fractional derivative of order a , with m — l<a<m, 
m e N, 

m 

D'^ fit) = D'^ g{t) ^ fit) = g{t) + Cj f^'' , (1-23) 

m 

D': fit) = D: git) ^ fit) = git) + J2 t^-' ■ (1-24) 

In these formulas the coefficients Cj are arbitrary constants. 

Incidentally, we note that (1.22) provides an instructive example to show how D°' 
is not right-inverse to J'* , since 

jaj^a^a-l = ^ ^ ^ ^ q ^ t > 0, . (1.25) 

For the two definitions we also note a difi'erence with respect to the formal limit 
as q; ^ (m- 1)+. From (1.13) and (1.17) we obtain respectively, 

a ^ (m - 1)+ =^ fit) D"^ J fit) = D"^'^ fit) ; (1.26) 

c^im-l)+ ^ D'^fit)^JD^fit) = D^-' fit) - /("^-i) (0+) . (1.27) 

We now consider the Laplace transform of the two fractional derivatives. For the 
standard fractional derivative the Laplace transform, assumed to exist, requires 
the knowledge of the (bounded) initial values of the fractional integral J"^~" and of 
its integer derivatives of order A; = 1, 2, m — 1 , as we learn from [2], [5], [10]. The 
corresponding rule reads, in our notation, 

m—l 

D'* /(t) ^ s'^ /(s) - ^ D'^ J^"^-'*^ /(0+) m-l<a<m. (1.28) 

fc=0 

The Caputo fractional derivative appears more suitable to be treated by the 
Laplace transform technique in that it requires the knowledge of the (bounded) ini- 
tial values of the function and of its integer derivatives of order A; = l,2,m — l,in 
analogy with the case when a = m . In fact, by using (1.10) and noting that 

m — l A. 

r fit) = J« J""-"" D'^ fit) = D"^ fit) = fit) - V Z^'^) (0+) - . (1.29) 

K ! 

fc=0 

we easily prove the following rule for the Laplace transform, 

m — l 

^* /(t) -«"/(«)- 5^/^'^ (0+)s"-'-^ m-l<a<m, (1.30) 

Indeed, the result (1.30), first stated by Caputo [19] by using the Fubini-Tonelli 
theorem, appears as the most "natural" generalization of the corresponding result 
well known for a = m . 
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We now show how both the definitions (1.13) and (1.17) for the fractional deriva- 
tive of f{t) can be derived, at least formally, by the convolution of $_Q,(t) with f{t) , 
in a sort of analogy with (1.8) for the fractional integral. For this purpose we need 
to recall from the treatise on generalized functions by Gel'fand and Shilov [16] that 
(with proper interpretation of the quotient as a limit if t = 0) 

^ — n— 1 

^-nit):=Y^ = S^''Ht), n = 0,l,... (1.31) 

where 5^'^\t) denotes the generalized derivative of order n of the Dirac delta distribu- 
tion. Here, we assume that the reader has some minimal knowledge concerning these 
distributions, sufficient for handling classical problems in physics and engineering. 

The equation (1.31) provides an interesting (not so well known) representation 
of 5^"^ (t) , which is useful in our following treatment of fractional derivatives. In 
fact, we note that the derivative of order n of a causal function f{t) can be obtained 
formally by the (generalized) convolution between and / , 

^f{t) = f^^\t) = ^_^{t)*f{t) = J^_f{T)S^^\t-T)dT, t>0, (1.32) 

based on the well known properties 




/(T)5^")(T-t)dT= (-1)"/("H^), S^^'^t - t) = S^^'^T - t) . (1.33) 



According to a usual convention, in (1.32-33) the limits of integration are extended 
to take into account for the possibility of impulse functions centred at the extremes. 
Then, a formal definition of the fractional derivative of order a could be 

The formal character is evident in that the kernel $_a(t) turns out to be not locally 
absolutely integrable and consequently the integral is in general divergent. In order 
to obtain a definition that is still valid for classical functions, we need to regularize the 
divergent integral in some way. For this purpose let us consider the integer m G N 
such that m — 1 < a < m and write —a — —m + {m — a) or —a — {m — a) — m. We 
then obtain 

[$_^(t) * ^m-a{t)] * f{t) = ^-m{t) * [^m-a{t) * f{t)] = D^J^-^f{t) , (1.34) 

or 

[^m-a{t) * ^-m{t)] * f{t) = $^-a(t) * [$_^(t) * f{t)] = J'— fit) . (1.35) 

As a consequence we derive two alternative definitions for the fractional derivative, 
corresponding to (1.13) and (1.17), respectively. The singular behaviour of ^-m{t) 
is refiected in the non-commutativity of convolution in these formulas. 
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1.4 Other Definitions and Notations 

Up to now we have considered the approach to fractional calculus usually re- 
ferred to Riemann and Liouville. However, while Riemann (1847) had generalized 
the integral Cauchy formula with starting point t = as reported in (1.1), originally 
Liouville (1832) had chosen t = — oo . In this case we define 

^-oo/W:= J-T r {t-Tr-'f{T)dT, ae]R+, (1.36) 

and consequently, form-l<Q!<m, meN, D"^^ f{t) := D"^ J-^"" f{t) , namely 

{rf"^ r 1 /■* /(r)rfT 1 
— r / 1 ^ — n ) m — 1 < a < m , 
dt^ [r{m-a)J_Ut-rr+^—\ ^^ ^^^ 

In this case, assuming f(t) to vanish as t — > — cxo along with its first m — 1 derivatives, 
we have the identity 

in contrast with (1.18). 

While for the fractional integral (1.2) a sufficient condition that the integral con- 
verge is that 

/(t) = 0(t'-^) , e>0, t^0+, (1.39) 
a sufficient condition that (1.36) converge is that 

/(t) = O (|t|-"-^) , e>0, t^-oo. (1.40) 

Integrable functions satisfying the properties (1.39) and (1.40) are sometimes referred 
to as functions of Riemann class and Liouville class, respectively [10]. For example 
power functions with 7 > — 1 and t > (and hence also constants) are of Riemann 
class, while \t\~^ with 5 > a > and t < and exp (ct) with c > are of Liouville 
class. For the above functions we obtain (as real versions of the formulas given in 
[10]) 

and 

J«^e^* = c-'^e'^*, D^o^e'^* = c'^e^*. (1.42) 
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Causal functions can be considered in the above integrals with the due care. In 
fact, in view of the possible jump discontinuities of the integrands at t = , in this 
case it is worthwhile to write 

{...)dT= [\...)dT. 



As an example we consider for < a < 1 the chain of identities 



1 /■* 

;i - a) Jo- 



fir) ^""_/(0+^ ■ 1 / 

a;) 



T{l-a) Jo- {t-r)"^ r(l-a)-" ' T{1 - a) Jo {t - t) 

f{0+)+D^f{t)=D'^f{t) 



(1.43) 



(1.45) 



T{l-a) 

where we have used (1.19) with m = 1 . 

In recent years it has become customary to use in place of (1.36) the Weyl frac- 
tional integral 

1 r°° 

WZfit):^TTr^ / (T-t)«-V(T)rfT, ae]R+, (1.44) 

based on a definition of Weyl (1917). For t > it is a sort of complementary integral 
with respect to the usual Riemann-Liouville integral (1.2). The relation between 
(1.36) and (1.44) can be readily obtained by noting that, see e.g. [10], 

f(t) = wzs /(^) = -f^ / + ^')""' /(-^') 

T{a) J_^ T{a) J^ 

= f{a) i ir'-tr-'fi-r')dT' = WZ9{t'), 

with g{t') = f{—t') , = — t . In the above passages we have carried out the changes 
of variable r — > r' = — r and t ^ t' = —t . 

For convenience of the reader, let us recall that exhaustive tables of Riemann- 
Liouville and Weyl fractional integrals are available in the second volume of the 
Bateman Project collection of Integral Transforms [16], in the chapter XIII devoted 
to fractional integrals. 

Last but not the least, let us consider the question of notation. The present 
authors oppose to the use of the notation for denoting the fractional integral, 

since it is misleading, even if it is used in distinguished treatises as [2], [10], [15]. It 
is well known that derivation and integration operators are not inverse to each other, 
even if their order is integer, and therefore such unification of symbols, present only in 
the framework of the fractional calculus, appears not justified. Furthermore, we have 
to keep in mind that for fractional order the derivative is yet an integral operator, so 
that, perhaps, it would be less disturbing to denote our as J~", than our J" as 
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1.5 The Law of Exponents 

In the ordinary calculus the properties of the operators of integration and differ- 
entiation with respect to the laws of commutation and additivity of their (integer) 
exponents are well known. Using our notation, the (trivial) laws 

where m, n = 0, 1, 2, . . . , can be referred to as the Law of Exponents for the operators 
of integration and differentiation of integer order, respectively. Of course, for any 
positive integer order, the operators and J"^ do not commute, see (1.11-12). 

In the fractional calculus the Law of Exponents is known to be generally true for 
the operators of fractional integration thanks to their semigroup property (1.3). In 
general, both the operators of fractional differentiation, and D'^ , do not satisfy 
either the semigroup property, or the (weaker) commutative property. To show how 
the Law of Exponents does not necessarily hold for the standard fractional derivative, 
we provide two simple examples (with power functions) for which 

(a) D^D^f{t)=D^D'^f{t)y^D"+^f{t), 

(b) D''D^g{t)^D'^D''g{t) = D^+^g{t). 

In the example (a) let us take f{t) = t"-*^/^ and a = P = 1/2. Then, using (1.15), 
we get D^/'^fit) = 0, D^/'^D^/'' f(t) = , but DV2+i/2y(^) = ^j(^) = -t-3/2/2 . 
In the example (b) let us take g{t) = t^^^ and a = 1/2, /? = 3/2. Then, again 
using (1.15), we get D^/^ g{t) = ^/7T/2, D^/^ g{t) = , but D^/'^D^/^ g{t) = 0, 
D^/^D^/^g{t) = -t3/2/4 and D^/^+^/^ g{t) = g{t) = -^3/2/4. 

Although modern mathematicians would seek the conditions to justify the Law of 
Exponents when the order of differentiation and integration are composed together, 
we resist the temptation to dive into the delicate details of the matter, but rather 
refer the interested reader to §IV.6 ("The Law of Exponents") in the book by Miller 
and Ross [10]. Let us, however, extract (in our notation, writing J'^ in place of D~°' 
for a > 0) three important cases, contained in their Theorem 3: If f{t) = t^ r]{t) or 
f{t) = t^lnt r]{t) , where A > — 1 and r]{t) = X^^o'^"^'^ having a positive radius R 
of convergence, then for < t < R , the following three formulas are valid: 

//>0 and {)<v<^i^D''J^f{t) = J^-''f{t), 
fi>0 and v> 11 ^ D'^J^ f{t)=D''-'' f{t), (1.48) 
0</x<A + l and > =^ D'' D^" f (t) = 0^"+" f (t) . 



At least in the case of f{t) without the factor Int, the proof of these formulas is 
straightforward. Use the definitions (1.2) and (1.13) of fractional integration and 
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differentiation, the semigroup property (1.3) of fractional integration, and apply the 
formulas (1.4) and (1.15) termwise to the infinite series you meet in the course of 
calculations. Of course, the condition that the function r]{t) be analytic can be 
considerably relaxed; it only need be "sufficiently" smooth." 

The lack of commutativity and the non-validity of the law of exponents has led 
to the notion of sequential fractional differentiation in which the order in which 
fractional differentiation operators D"^ , ^ . . . ^ are concatenated is crucial. 
For this and the related field of fractional differential equations we refer again to 
Miller and Ross [10]. Furthermore, Podlubny [29] has also given formulas for the 
Laplace transforms of sequential fractional derivatives. 

In order to give an impression on the strange effects to be expected in use of 
sequential fractional derivatives we consider for a function f{t) continuous for t > 
and for positive numbers a and P with a + P = 1 the three problems (a), (b), (c) with 
the respective general solutions u,v ,w in the set of locally integrable functions, 

(a) D^D^ u{t) = fit) u{t) = J fit) +ai+a2 1^'^ , 

(b) Df^D"" v{t) = fit) => vit) = J fit) + bi + b2 f^-^ , (1.49) 
ic)Duit) = fit) wit) = Jfit) + c, 

where ai , a2 , 6i , 62 , c are arbitrary constants. Whereas the result for (c) is obvious, 
in order to obtain the final results for (a) [or (b)] we need to apply first the operator 
J" [or J^] and then the operator [or J"]. The additional terms must be taken 
into account because t'^~^ = 0, J^~'^ t^~^ = r(7) , 7 = 0;,/?. We observe that, 
whereas the general solution of (c) contains one arbitrary constant, that of (a) and 
likewise of (6) contains two arbitrary constants, even though a + (3 = 1 .In case (3 
the singular behaviour of uit) at t = O"*" is distinct from that of vit) . 

From above we can conclude in rough words: sufficiently fine sequentialization 
increases the number of free constants in the general solution of a fractional differ- 
ential equation, hence the number of conditions that must be imposed to make the 
solution unique. For an example see Bagley's treatment of a composite fractional os- 
cillation equation [30]; there the highest order of derivative is 2, but four conditions 
are required to achieve uniqueness. 

In the present lectures we shall avoid the above troubles since we shall consider 
only differential equations containing single fractional derivatives. Furthermore we 
shall adopt the Caputo fractional derivative in order to meet the usual physical re- 
quirements for which the initial conditions are expressed in terms of a given number 
of bounded values assumed by the field variable and its derivatives of integer order, 
see (1.24) and (1.30). 
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2. FRACTIONAL INTEGRAL EQUATIONS 

In this section we shall consider the most simple integral equations of fractional 
order, namely the Abel integral equations of the first and the second kind. The former 
investigations on such equations are due to Abel (1823-26), after whom they are 
named, for the first kind, and to Hille and Tamarkin (1930) for the second kind. The 
interested reader is referred to [5], [21] and [31-33] for historical notes and detailed 
analysis with applications. Here we limit ourselves to put some emphasis on the 
method of the Laplace transforms, that makes easier and more comprehensible the 
treatment of such fractional integral equations, and provide some applications. 

2.1 Abel integral equation of the first kind 

Let us consider the Abel integral equation of the first kind 
1 /"* u(t) 

where f{t) is a given function. We easily recognize that this equation can be expressed 
in terms of a fractional integral, i. e. 

ru{t)^f{t), (2.2) 

and consequently solved in terms of a fractional derivative, according to 

u{t)=D^f{t). (2.3) 

To this end we need to recall the definition (1.2) and the property (1.14) J"' = I. 

Let us now solve (2.1) using the Laplace transform. Noting from (1.7-8) and 
(1.10) that J"u(t) = $a(^) * u{t) ^ u{s)/s'^ , we then obtain 

^ = f{s)^u{s) = s-f{s). (2.4) 

Now we can choose two different ways to get the inverse Laplace transform from 
(2.4), according to the standard rules. Writing (2.4) as 



u{s) = s 

we obtain 



M 

gl-a 



(2.4a) 



u(t) = — r 4- f /^^\ dr. (2.5a) 
^ ^ r(l - a) dt Jo {t - t)« ^ ^ 



On the other hand, writing (2.4) as 



1 , /(O^) 



s 

we obtain 



u(t) = / , ^'^^l dr + /(0+) " , . (2.56) 

r(l - a) Jo (t-r)^ •'^ ^ r(l -a) ^ ^ 
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Thus, the solutions (2.5a) and (2.5b) are expressed in terms of the fractional deriva- 
tives D"' and , respectively, according to (1.13), (1.17) and (1.19) with m = 1 . 

The way b) requires that f{t) be differentiable with >C-transformable derivative; 
consequently < 1/(0+) | < oo . Then it turns out from (2.5b) that tt(0"'") can be 
infinite if /(O"*") 7^ , being u(t) = 0(t~") , as t ^ 0+ . The way a) requires weaker 
conditions in that the integral at the R.H.S. of (2.5a) must vanish as t — 0"*"; conse- 
quently /(0+) could be infinite but with f{t) = 0{t~'') , 0<z/<l-a ast^O+. 
To this end keep in mind that $i-q * ^i-v — ^2-ci-v Then it turns out from (2.5a) 
that m(0+) can be infinite if /(0+) is infinite, being u{t) = 0(t-^°=+'')) , as t ^ 0+ . 

Finally, let us remark that we can analogously treat the case of equation (2.1) 
with < CK < 1 replaced hya>0.1fm — l<a<m with m e N , then again we 
have (2.2), now with f{t) given by the formula (1.13) which can also be obtained 
by the Laplace transform method. 

2.2 Abel integral equation of the second kind 

Let us now consider the Abel equation of the second kind 




(2.7) 



(2.6) 



u 



{t) = (1 + AJ«)-^ fit) = 1 + 5](-A)- J"" fit) . 




(2.8) 



Noting by (1.7-8) that 



J""/(t) = $an(t) * fit) 



ctn — 1 

+ 



* fit) 



r(an) 



the formal solution reads 




(2.9) 



Recalling from the Appendix the definition of the function. 




t>0, a>0, AgC, 



(2.10) 



where denotes the Mittag-Leffler function of order a , we note that 




(2.11) 
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Finally, the solution reads 

u{t) = fit) + e'^{t;X) * fit). (2.12) 

Of course the above formal proof can be made rigorous. Simply observe that 
because of the rapid growth of the gamma function the infinite series in (2.9) and 
(2.11) are uniformly convergent in every bounded interval of the variable t so that 
term-wise integrations and differentiations are allowed. However, we prefer to use 
the alternative technique of Laplace transforms, which will allow us to obtain the 
solution in different forms, including the result (2.12). 

Applying the Laplace transform to (2.6) we obtain 



1+4' 



u{s) = f{s) 



u{s) 



(2.13) 



s« + A 

Now, let us proceed to obtain the inverse Laplace transform of (2.13) using the 
following Laplace transform pair (see Appendix) 

„a-l 



(2.14) 



s« + A 

As for the Abel equation of the first kind, we can choose two different ways to get 
the inverse Laplace transforms from (2.13), according to the standard rules. Writing 
(2.13) as 



u{s) 



+ A 



we obtain 

If we write (2.13) as 

u{s) 

we obtain 

uit) 



d /■* 

«W = ^ JJ{t-T)ea{T;X)dT. 



+ A 

ft 



[^/»-/(o+)] + /(o+) 



+ A 



(2.13a) 
(2.15a) 

(2.136) 



/ fit - t) cair; A) dr + /(0+) e^it; A) . (2.156) 
Jo 

We also note that, e^it; X) being a function differentiable with respect to t with 
60,(0+; A) = £'0,(0+) = 1 , there exists another possibility to re- write (2.13), namely 

„a-l 



u{s) = 



s« + A 



- 1 



Then we obtain 



u{t)= ff{t-T)e'^{T-X)dT + f{i) 
Jo 



(2.13c) 



(2.15c) 



in agreement with (2.12). We see that the way h) is more restrictive than the ways a) 
and c) since it requires that fit) be differentiable with >C-transformable derivative. 
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2.3 Some applications of Abel integral equations 

It is well known that Niels Henrik Abel was led to his famous equation by the 
mechanical problem of the tautochrone, that is by the problem of determining the 
shape of a curve in the vertical plane such that the time required for a particle to 
slide down the curve to its lowest point is equal to a given function of its initial height 
(which is considered as a variable in an interval [0, ff]). After appropriate changes 
of variables he obtained his famous integral equation of first kind with a = 1/2 . He 
did, however, solve the general case < a < 1 . See Tamarkin's translation of and 
comments to Abel's short paper . As a special case Abel discussed the problem of 
the isochrone, in which it is required that the time of sliding down is independent of 
the initial height. Already in his earlier publication he recognized the solution as 
derivative of non- integer order. 

We point out that integral equations of Abel type, including the simplest (2.1) and 
(2.6), have found so many applications in diverse fields that it is almost impossible 
to provide an exhaustive list of them. 

Abel integral equations occur in many situations where physical measurements 
are to be evaluated. In many of these the independent variable is the radius of a 
circle or a sphere and only after a change of variables the integral operator has the 
form J" , usually with a = 1/2 , and the equation is of first kind. Applications are, 
e.g. , in evaluation of spectroscopic measurements of cylindrical gas discharges, the 
study of the solar or a planetary atmosphere, the investigation of star densities in 
a globular cluster, the inversion of travel times of seismic waves for determination 
of terrestrial sub-surface structure, spherical stereology. Descriptions and analysis of 
several problems of this kind can be found in the books by Gorenflo and Vessella [21] 
and by Craig and Brown [31], see also [32]. Equations of first and of second kind, 
depending on the arrangement of the measurements, arise in spherical stereology. 
See [33] where an analysis of the basic problems and many references to previous 
literature are given. 



Abel, N. H.: Solution of a mechanical problem. Translated from the German. In: D. E. Smith, 
editor: A Source Book in Mathematics, pp. 656-662. Dover Publications, New York, 1959. 

^■^ Abel. N. H.: Auflocsung cincr mcchanischen Aufgabe. Journal fiir die reine und angewandte 

Mathcmatik (CrcUc), Vol. I (1826), pp. 153-157. 

^) Abel, N. H.: Solution de quelques problemes a I'aide d'integrales definies. Translated from the 
Norwegian original, published in Magazin for Naturvidenskaberne. Aargang 1, Bind 2, Christiana 
1823. French Translation in Oeuvres Completes, Vol I, pp. 11-18. Nouvelle edition par L. Sylow 
et S. Lie, 1881. 
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Another field in which Abel integral equations or integral equations with more 
general weakly singular kernels are important is that of inverse boundary value prob- 
lems in partial difi^erential equations, in particular parabolic ones in which naturally 
the independent variable has the meaning of time. We are going to describe in detail 
the occurrence of Abel integral equations of first and of second kind in the problem 
of heating (or cooling) of a semi-infinite rod by infiux (or efflux) of heat across the 
boundary into (or from) its interior. Consider the equation of heat flow 

Ut-Uxx^O, u^u{x,t), (2.16) 

in the semi-infinite intervals < a; < cxd and < t < oo of space and time, re- 
spectively. In this dimensionless equation u = u{x,t) means temperature. Assume 
vanishing initial temperature, i.e. u{x^ 0) = for < a; < oo and given infiux across 
the boundary a; = from a; < to a; > , 

-Ux{Q,t)=p{t). (2.17) 

Then, under appropriate regularity conditions, ■u(a;,t) is given by the formula, see 
e.g. [34], 



^a;,t) = ^ r4^e-^V[4(*-T)]^^^ a; > , t > . (2.18) 



We turn our special interest to the interior boundary temperature (j){t) := u{0~^,t) , 
t > , which by (2.18) is represented as 

J- f J^dT = J^^^p{t)=m. ^>0. (2.19) 

We recognize (2.19) as an Abel integral equation of first kind for determination of an 
unknown influx p{t) if the interior boundary temperature (p{t) is given by measure- 
ments, or intended to be achieved by controlling the infiux. Its solution is given by 
formula (1-13) with m — 1 , a — 1/2 , as 

It may be illuminating to consider the following special cases, 

(i) cf>{t) = t ^ p{t) = lv^t, 

(2.21) 

(n) 4,{t) = 1 =^ p{t) = ^ 



'TTt 

where we have used formula (1.15). So, for linear increase of interior boundary 
temperature the required infiux is continuous and increasing from towards oo (with 
unbounded derivative at t = 0+), whereas for instantaneous jump- like increase from 
to 1 the required infiux decreases from oo at t = 0+ to as t — > oo . 
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We now modify our problem to obtain an Abel integral equation of second kind. 
Assume that the rod a; > is bordered at a; = by a bath of liquid in a; < with 
controlled exterior boundary temperature u{0~ ,t) := '4>{t) . 

Assuming Newton's radiation law we have an influx of heat from 0~ to O"*" pro- 
portional to the difference of exterior and interior temperature, 

p{t) = X [iP(t) - (f)(t)] , A > . (2.22) 

Inserting (2.22) into (2.19) we obtain 

VTT Jo y/t-T 

namely, in operator notation, 

(l + X (j){t) = X ^(^) . (2.23) 

If we now assume the exterior boundary temperature '4>{t) as given and the evolution 
in time of the interior boundary temperature (^{t) as unknown, then (2.23) is an Abel 
integral equation of second kind for determination of (f>{t) . 

With a = 1/2 the equation (2.23) is of the form (2.7), and by (2.8) its solution is 

4>{t) = A (l + A " ^(^) {-xr+^ j(-+i)/2 ^(^) . (2.24) 

m=0 

Let us investigate the very special case of constant exterior boundary temperature 

^{t) = l. (2.25) 

Then, by (1.4) with 7 = 0, 

f{m+l)/2 

^^^^-r[(m + l)/2 + l]' 

hence 

°° j.(m+l)/2 oo j-n/2 

,(0 = - 5:(-A)"- = 1 - E(-A)" , 

m—0 '■^ ^' ^ n=0 ^ ' ^ 

SO that 

cf>{t) = 1 - (-At^/^) = 1 - ei/2(t; A) . (2.26) 

Observe that 0(t) is creep function, increasing strictly monotonically from towards 
1 as t runs from to oo . 

For more or less distinct treatments of this problem of "Newtonian heating" the 
reader may consult [21], and [35-37]. In [37] a formulation in terms of fractional dif- 
ferential equations is derived and, furthermore, the analogous problem of "Newtonian 
cooling" is discussed. 
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3. FRACTIONAL DIFFERENTIAL EQUATIONS: 1-st PART 

We now analyse the most simple differential equations of fractional order. For this 
purpose, following our recent works [37-42] , we choose the examples which, by means 
of fractional derivatives, generalize the well-known ordinary differential equations 
related to relaxation and oscillation phenomena. In this section wc treat the simplest 
types, which we refer to as the simple fractional relaxation and oscillation equations. 
In the next section we shall consider the types, somewhat more cumbersome, which 
we refer to as the composite fractional relaxation and oscillation equations. 

3. 1 The simple fractional relaxation and oscillation equations 

The classical phenomena of relaxation and oscillations in their simplest form are 
known to be governed by linear ordinary differential equations, of order one and two 
respectively, that hereafter we recall with the corresponding solutions. Let us denote 
by -u = u{t) the field variable and by q{t) a given continuous function, with t > . 

The relaxation differential equation reads 

u'{t) = -u{t)+q{t), (3.1) 

whose solution, under the initial condition tt(0"'") = co , is 

j-t 

u{t) = coe~^+ q{t-T)e~^dT. (3.2) 
Jo 

The oscillation differential equation reads 

u"{t) = -u{t) + q{t) , (3.3) 
whose solution, under the initial conditions u{0~^) — cq and u'{0'^) = ci , is 

u{t) = Co cost -|- ci sint -|- q{t — r) sin r dr. (3-4) 

Jo 

From the point of view of the fractional calculus a natural generalization of equa- 
tions (3.1) and (3.3) is obtained by replacing the ordinary derivative with a fractional 
one of order a . In order to preserve the type of initial conditions required in the clas- 
sical phenomena, we agree to replace the first and second derivative in (3.1) and (3.3) 
with a Caputo fractional derivative of order a with < a < 1 and 1 < a < 2 , re- 
spectively. We agree to refer to the corresponding equations as the simple fractional 
relaxation equation and the simple fractional oscillation equation. 
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Generally speaking, we consider the following differential equation of fractional 
order a > , 

"(^)-E^^^'^(O^) =-^(^)+«(^)' ^>0- (3.5) 
fc=o ■ / 

Here m is a positive integer uniquely defined hym— l<a<m, which provides the 
number of the prescribed initial values u^'^\0'^) = Ck , A; = 0, 1, 2, . . . , m — 1 . Implicit 
in the form of (3.5) is our desire to obtain solutions u{t) for which the u^'^^t) are 
continuous for t > 0, k = 0, l,...,m — 1. In particular, the cases of fractional 
relaxation and fractional oscillation are obtained for m = 1 and m = 2 , respectively 
We note that when a is the integer m the equation (3.5) reduces to an ordinary 
differential equation whose solution can be expressed in terms of m linearly inde- 
pendent solutions of the homogeneous equation and of one particular solution of the 
inhomogeneous equation. We summarize the well-known result as follows 

m-l „i 

^(*) = '^CkUk{t)+ / q{t - t) us{t) dr . (3.6) 
fe=o 

Ukit) = j'' uoit) , ui^\0+) = Skh, h,k = 0,l,...,m-l, (3.7) 

U5it) = -u'o{t). (3.8) 

Thus, the m functions Uk{t) represent the fundamental solutions of the differential 
equation of order m , namely those linearly independent solutions of the homoge- 
neous equation which satisfy the initial conditions in (3.7). The function us{t) , 
with which the free term q{t) appears convoluted, represents the so called impulse- 
response solution, namely the particular solution of the inhomogeneous equation with 
all Cfc = , k — 0,1, . . . ,m — 1 , and with q(t) = S{t) . In the cases of ordinary re- 
laxation and oscillation we recognize that uo{t) = e~* = us{t) and UQ{t) — cos t, 
ui{t) — Juo{t) = sin t = cos (t — 7r/2) = us{t) , respectively. 
Remark 1 : The more general equation 

/ "^-1 fk \ 

i;"U(t)-5^-«('=)(0+) =-p"^t) + g(t), p>0, t>0, (3.50 
\ fc=o ■ / 

can be reduced to (3.5) by a change of scale t ^ t/p . We prefer, for ease of notation, 
to discuss the " dimensionless" form (3.5). 

Let us now solve (3.5) by the method of Laplace transforms. For this purpose 
we can use directly the Caputo formula (1.30) or, alternatively, reduce (3.5) with 
the prescribed initial conditions as an equivalent (fractional) integral equation and 
then treat the integral equation by the Laplace transform method. Here we prefer 
to follow the second way. Then, applying the operator J'* to both sides of (3.5) we 
obtain 
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m—1 1^ 

u{t) = E - + ^(^) • (3.9) 



A;=0 

The application of the Laplace transform yields 

fc=0 



hence 

^(^) = E + ^~(^) • (3-10) 

Introducing the Mittag-LefHer type functions 



e^{t) = e«(t; 1) := E^{-t^) - ^— - , (3.11) 

«fc(t) := J'=ea(t) - A; = 0,l,...,m-1, (3.12) 

we find, from inversion of the Laplace transforms in (3.10), 

m-l „t 

= E'^'=^'=(^)~ / <7(^-'r)«o('^)c^T. (3.13) 

fc=0 

For finding the last term in the R.H.S. of (3.13), we have used the well-known rule 
for the Laplace transform of the derivative noting that uo(0"'") = 6^(0+) = 1 , and 

1 / s"-^ 

+ 1 " 

The formula (3.13) encompasses the solutions (3.2) and (3.4) found for a = 1 , 2 , 
respectively. 

When a is not integer, namely for m — 1 < a < m , we note that m — 1 represents 
the integer part of a (usually denoted by [a] ) and m the number of initial conditions 
necessary and sufficient to ensure the uniqueness of the solution u{t). Thus the m 
functions Uk{t) = J^eo.{t) with A; = 0, 1, ... , m—1 represent those particular solutions 
of the homogeneous equation which satisfy the initial conditions 

(0+) = 5uh, /i, = 0, 1, . . . , m - 1 , (3.15) 



(^^-1) ^ -u',{t) = -e'^{t) . (3.14) 



and therefore they represent the fundamental solutions of the fractional equation 
(3.5), in analogy with the case a = m. Furthermore, the function us{t) = —e'^{t) 
represents the impulse-response solution. Hereafter, we are going to compute and 
exhibit the fundamental solutions and the impulse-response solution for the cases (a) 
< ct < 1 and (b) 1 < a < 2 , pointing out the comparison with the corresponding 
solutions obtained when a = 1 and a = 2 . 



244 Fractional Calculus: Integral and Differential Equations of Fractional Order 



Remark 2 : The reader is invited to verify that the solution (3.13) has continuous 
derivatives u^''\t) for k = 0, 1, 2, . . . , m — 1 , which fulfill the m initial conditions 
ix(^)(0+) = Cfc . In fact, looking back at (3.9), one must recognize the smoothing 
power of the operator J" . However, the so called impulse-response solution of our 
equation (3.5), us{t) , is expected to be not so regular like the ordinary solution (3.13). 
In fact, from (3.10) and (3.13-14), one obtains 

usit) = -u'oit) ^ , (3.16) 

and therefore, using the limiting theorem for Laplace transforms, one can recognize 
that, being m — 1 < a < m , 

u'i'\0+) = 0, h = 0,l,...,m-2; u^p-^\o+) ^ oo . (3.17) 

We now prefer to derive the relevant properties of the basic functions Ca (t) directly 
from their representation as a Laplace inverse integral 

in detail for < o: < 2 , without detouring on the general theory of Mittag-Leffler 
functions in the complex plane. In (3.18) Br denotes the Bromwich path, i.e. a line 
Re {s} = a with a value cr > 1 , and Im{s} running from — oo to +00 . 
For transparency reasons, we separately discuss the cases 

(a) < a < 1 and (b) 1 < a < 2 , 

recalling that in the limiting cases a = 1 , 2 , we know ea{t) as elementary function, 
namely ei{t) = and 62 (t) = cos t. For a not integer the power function s" is 

uniquely defined as — e**^"^^* , with — tt < args < tt , that is in the complex 
s-plane cut along the negative real axis. 

The essential step consists in decomposing ea{t) into two parts according to 
Ga{t) — fa{t) + Qait) 7 as indicated below. In case (a) the function fa{t) , in case 
(b) the function —fait) is completely monotone; in both cases faif) tends to zero 
as t tends to infinity, from above in case (a), from below in case (b). The other 
part, ga{t) , is identically vanishing in case (a), but of oscillatory character with 
exponentially decreasing amplitude in case (b). 

In order to obtain the desired decomposition of we bend the Bromwich path 
of integration Br into the equivalent Hankel path Ha{l~^), a loop which starts from 
—00 along the lower side of the negative real axis, encircles the circular disc |s| = 1 
in the positive sense and ends at — oo along the upper side of the negative real axis. 
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One obtains 



eait) = fait) + gait) , t>0, 



with 



fait) 



.St 



Ha(e) 



+ 1 



ds 



(3.19) 



(3.20) 



where now the Hankel path Haie) denotes a loop constituted by a small circle \s\ = e 
with e ^ and by the two borders of the cut negative real axis, and 



9ait) :=5^e^''»*i2es 



+ 1 



1 



a ^ 



(3.21) 



where s'f^ are the relevant poles of s'*~-^/(s" + 1). In fact the poles turn out to be 
Sh = exp [z(2/j, + 1)77/0;] with unitary modulus; they are all simple but relevant are 
only those situated in the main Riemann sheet, i.e. the poles s'f^ with argument such 
that — TT < arg < tt . 

If < a < 1 , there are no such poles, since for all integers h we have | arg Sh\ = 
\2h + 1\ n / a > n ; as a consequence, 



g^t) = , hence Cait) = fait) , if < a < 1 



(3.22) 



If 1 < q; < 2 , then there exist precisely two relevant poles, namely Sq = exp(z7r/Q;) 
and s'_i = exp(— ztt/q:) = sq' , which are located in the left half plane. Then one 
obtains 

t cos (tt/q!) 



9ait) 



COS 



a 



t sin 



-) 

a/ 



if 1< a < 2 . 



(3.23) 



We note that this function exhibits oscillations with circular frequency a; (a) = 
sin (tt/q;) and with an exponentially decaying amplitude with rate A (a) = | cos (tt/q;) | . 

Remark 3 : One easily recognizes that (3.23) is valid also for 2 < a < 3 . In the 
classical case a = 2 the two poles are purely imaginary (coinciding with ±z) so that 
we recover the sinusoidal behaviour with unitary frequency. In the case 2 < ct < 3 , 
however, the two poles are located in the right half plane, so providing amplified 
oscillations. This instability, which is common to the case a = 3 , is the reason why 
we limit ourselves to consider a in the range < ct < 2 . 

It is now an exercise in complex analysis to show that the contribution from the 
Hankel path Haie) as e — > is provided by 



/•oo 

fait):= / e-^^Ka(r)dr, 
Jo 



with 



Kair) = Im 

TT 



,a-l 



+ 1 



.a-1 



sin (an) 



TT r^" + 2 r'* cos (cktt) + 1 



(3.24) 



(3.25) 
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This function Ka{r) vanishes identically if a is an integer, it is positive for all r 
if < a < 1 , negative for all r if 1 < a < 2 . In fact in (3.25) the denominator is, 
for a not integer, always positive being > (r'* — 1)^ > . Hence fa{t) has the afore- 
mentioned monotonicity properties, decreasing towards zero in case (a), increasing 
towards zero in case (b). We also note that, in order to satisfy the initial condition 
e„(0+) = 1, we find Ka{r) dr = lifO<a<l, Ka{r) dr = 1 - 2/a if 
1 < q: < 2 . In Fig. 1 we show the plots of the spectral functions Kq, (r) for some 
values of a in the intervals (a) < a < 1 , (b) 1 < o; < 2 . 
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Fig. la - Plots of the basic spectral function Kq, (r) for < a < 1 
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Fig. lb - Plots of the hasic spectral function —Ka{r) for 1 < a < 2 
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In addition to the basic fundamental solutions, uo{t) = ea{t) we need to compute 
the impulse-response solutions us{t) = —D^ ea{t) for cases (a) and (b) and, only in 
case (b), the second fundamental solution ui{t) = ea{t) . For this purpose we note 
that in general it turns out that 



with 



J" fait) 



-i-k 



sin (an) 



(3.26) 



(3.27) 



(3.27) 



r^" + 2 r" cos (cutt) + 1 ' 

where i^a(r) = Kafi{r) , and 

J^Oait) = - e^ (^/°^) cos \t sin ( -) - k- . 

a I \a/ a. 

This can be done in direct analogy to the computation of the functions ea{t), the 
Laplace transform of J^eait) being given by (3.12). For the impulse-response solution 
we note that the effect of the differential operator is the same as that of the virtual 
operator . 

In conclusion we can resume the solutions for the fractional relaxation and oscil- 
lation equations as follows: 

(a) < a < 1 , 

u{t) = coUo{t) + q{t - t) us{t) dr , 
Jo 



(3.28a) 



where 



Uq 



{t)= / 

^0 



e-^*K«,o(r)rfr, 



usit) = - / 

^0 

with uo{0+) = 1 , U5{0+) = oo ; 
(b) 1 < a < 2 , 



-rt 



(3.29a) 



Ka-i{r)dr, 



where 



u{t) = coUo{t) + ciui{t) + / q{t - t) us{t) dr , 

Jo 

' uo{t) = He K^,o{r) dr + - (^/«) 
Jo « 

u,{t)^ He 1 (r) dr+-et i^/a) 

Jo ' oi 

usit) = - / 



(3.286) 



t sin 
t sin 



TT 
CK/ J 
TT 

a 



a 



(3.296) 



e~" Ka-iir)dr- - e 



— — pt cos (tt/q;) 



+ 



TT 



a 



with «o(0+) = 1, ^0(0+) = 0, ^1(0+) = 0, u[iO+) = 1, «5(0+) = 0, ^^(0+) = +00. 
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Fig. 2b - Plots of the basic fundamental solution uo{t) = ea{t) for 1 < a <2 
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In Fig. 2 we quote the plots of the basic fundamental solution for the following 
cases : (a) a = 0.25, 0.50, 0.75, 1, and (b) a = 1.25, 1.50, 1.75, 2, obtained 
from the first formula in (3.29a) and (3.29b), respectively. We have verified that 
our present results confirm those obtained by Blank [43] by a numerical treatment 
and those obtained by Mainardi [39] by an analytical treatment, valid when a is a 
rational number, see §A2 of the Appendix. Of particular interest is the case a = 1/2 
where we recover a well-known formula of the Laplace transform theory, see (A. 34), 

ei/2(t) := E,/2{-Vi) = e^erfc(v^) - ^i/, ^j/, ^ , (3.30) 

where erfc denotes the complementary error function. 

We now desire to point out that in both the cases (a) and (b) (in which a is just not 
integer) i.e. for fractional relaxation and fractional oscillation, all the fundamental 
and impulse-response solutions exhibit an algebraic decay as t — > oo , as discussed 
below. Let us start with the asymptotic behaviour of uo{t) . To this purpose we first 
derive an asymptotic series for the function fa{t), valid for t ^ oo . Using the identity 



s« -Fl ' ' -M 

in formula (3.20) and the Hankel representation of the reciprocal Gamma function, 
we (formally) obtain the asymptotic expansion (for a non integer) 

N 

m ^ r(l -na) ^^ , as t ^ oc . (3.31) 

n=l ^ ^ 

The validity of this asymptotic expansion can be established rigorously using the 
(generalized) Watson lemma, see [44]. We also can start from the spectral represen- 
tation (3.24-25) and expand the spectral function for small r. Then the (ordinary) 
Watson lemma yields (3.31). We note that this asymptotic expansion coincides with 
that for Uo{t) — ea{t), having assumed 0<q:<2(q:7!^1). In fact the contribution of 
ga{t) is identically zero if < ct < 1 and exponentially small ast— >ooifl<Q;<2. 

The asymptotic expansions of the solutions ui{t) and us{t) are obtained from 
(3.31) integrating or differentiating term by term with respect to t . In particular, 
taking the leading term in (3.31), we obtain the asymptotic representations 

y^ — n: ^1 — a ^ — a — 1 

"°^'^^r(r^' "^^'^~f(^' ""'^'^^-fi^)' ^^-^'^ 

that point out the algebraic decay of the fundamental and impulse-response solutions. 

In Fig. 3 we show some plots of the basic fundamental solution UQ{t) = ea{t) 
for a = 1.25 , 1.50 , 1.75 . Here the algebraic decay of the fractional oscillation can 
be recognized and compared with the two contributions provided by /q, (monotonic 
behaviour ) and ga{t) (exponentially damped oscillation). 
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Fig. 3a - Decay of the basic fundamental solution UQ{t) = ea{t) for a = 1.25 



0.05 



-0.05 









0=1.50 
















t 


5 ' / 




) ^ ' 


I-? 


-.'' /■■ ' y** 









Fig. 3b - Decay of the basic fundamental solution uo{t) = ea{t) for a = 1.50 
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Fig. 3c - Decay of the basic fundamental solution uo{t) — ea{t) for a — 1.75 
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3. 2 The zeros of the solutions of the fractional oscillation equation 

Now we find it interesting to carry out some investigations about the zeros of the 
basic fundamental solution uo{t) = ea{t) in the case (b) of fractional oscillations. 
For the second fundamental solution and the impulse-response solution the analysis 
of the zeros can be easily carried out analogously. 

Recalling the first equation in (3.29b), the required zeros of ea{t) are the solutions 
of the equation 

ea{t) = fait) + - e^ (^/o^) COS sin = . (3.33) 

a I \a/i 

We first note that the function eo.{t) exhibits an odd number of zeros, in that 
ea(0) = 1 , and, for sufficiently large t, ea{t) turns out to be permanently negative, 
as shown in (3.32) by the sign of r(l — a) . The smallest zero lies in the first positivity 
interval of cos [t sin (tt/q;)] , hence in the interval < t < 7r/[2 sin (tt/q;)] ; all other 
zeros can only lie in the succeeding positivity intervals of cos [t sin (tt/q;)] , in each of 
these two zeros are present as long as 

2gtcos(7r/a) > _ (3 34) 

a 

When t is sufficiently large the zeros are expected to be found approximately from 
the equation 

2 cos (tt/q;) ^ ^ (3.35) 

a \T{1 — a)\ 

obtained from (3.33) by ignoring the oscillation factor of ga{t) [see (3.23)] and taking 
the first term in the asymptotic expansion of fa{t) [see (3.31-32)]. As we have shown 
in [40], such approximate equation turns out to be useful when o; — > l""" and o; — > 2~ . 

For CK — > 1+ , only one zero is present, which is expected to be very far from the 
origin in view of the large period of the function cos [t sin (tt/q;)] . In fact, since there 
is no zero for a = 1, and by increasing a more and more zeros arise, we are sure that 
only one zero exists for a sufficiently close to 1. Putting a = 1 + e the asymptotic 
position T* of this zero can be found from the relation (3.35) in the limit e ^ 0"*" . 
Assuming in this limit the first-order approximation, we get 

~ log 0^ , (3.36) 

which shows that tends to infinity slower than 1/e , as e — > . For details see [40] . 
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For CK — > 2~, there is an increasing number of zeros up to infinity since 62 (t) = cos t 
has infinitely many zeros [in t* = (n + l/2)7r , n= 0, 1, . . .]. Putting now a = 2 — 5 
the asymptotic position T* for the largest zero can be found again from (3.35) in the 
limit 5 — > 0+ . Assuming in this limit the first-order approximation, we get 

T.~iilog(i). (3.37) 

For details see [40]. Now, for 5—^0"'" the length of the positivity intervals of goi{t) 
tends to n and, as long as t < , there are two zeros in each positivity interval. 
Hence, in the limit 5—^0"'", there is in average one zero per interval of length tt , so 
we expect that A^* ~ T* /tt . 

Remark 4 : For the above considerations we got inspiration from an interesting paper 
by Wiman [45] who at the beginning of our century, after having treated the Mittag- 
Lefiier function in the complex plane, considered the position of the zeros of the 
function on the negative real axis (without providing any detail). Our expressions 
of Th, are in disagreement with those by Wiman for numerical factors; however, the 
results of our numerical studies carried out in [40] confirm and illustrate the validity 
of our analysis. 

Here, we find it interesting to analyse the phenomenon of the transition of the 
(odd) number of zeros as 1.4 < ct < 1.8 . For this purpose, in Table I we report the 
intervals of amplitude Aa = 0.01 where these transitions occur, and the location 
T* (evaluated within a relative error of 0.1%) of the largest zeros found at the two 
extreme values of the above intervals. We recognize that the transition from 1 to 3 
zeros occurs as 1.40 < a < 1.41, that one from 3 to 5 zeros occurs as 1.56 < o; < 1.57, 
and so on. The last transition in the considered range of a is from 15 to 17 zeros, 
and it just occurs as 1.79 < a < 1.80 . 





a 




1^3 


1.40^ 1.41 


1.730^5.726 


3 ^ 5 


1.56 ^ 1.57 


8.366 ^ 13. 18 


5^7 


1.64^ 1.65 


14.61 ^ 20.00 


7^9 


1.69^ 1.70 


20.80 ^ 26.33 


9^ 11 


1.72^ 1.73 


27.03 ^ 32.83 


11^13 


1.75 ^ 1.76 


33.11 ^38.81 


13^ 15 


1.78^ 1.79 


39.49 -^ 45.51 


15h- 17 


1.79 1.80 


45.51 ^51.46 



Table I 

= number of zeros, a = fractional order, location of the largest zero. 
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4. FRACTIONAL DIFFERENTIAL EQUATIONS: 2-nd PART 

In this section we shall consider the following fractional differential equations for 
t > , equipped with the necessary initial conditions, 

+ + m(0+) = co, 0<a<l, (4.1) 

^+a^+^(t)=g(^), v{0+)^co, v'{0+)^ci, 0<a<2, (4.2) 

where o is a positive constant. The unknown functions u{t) and v{t) (the field vari- 
ables) are required to be sufficiently well behaved to be treated with their derivatives 
u'{t) and v'{t) , v"{t) by the technique of Laplace transform. The given function q{t) 
is supposed to be continuous. In the above equations the fractional derivative of 
order a is assumed to be provided by the operator the Caputo derivative, see 
(1.17), in agreement with our choice followed in the previous section. Note that in 
(4.2) we must distinguish the cases (a) < a < 1 , (b) 1 < a < 2 and a = 1 . 

The equations (4.1) and(4.2) will be referred to as the composite fractional relax- 
ation equation and the composite fractional oscillation equation, respectively, to be 
distinguished from the corresponding simple fractional equations treated in §3. 

The fractional differential equation in (4.1) with a = 1/2 corresponds to the 
Basset problem, a classical problem in fluid dynamics concerning the unsteady motion 
of a particle accelerating in a viscous fluid under the action of the gravity, see [24] . 

The fractional differential equation in (4.2) with < ct < 2 models an oscillation 
process with fractional damping term. It was formerly treated by Caputo [19], who 
provided a preliminary analysis by the Laplace transform. The special cases a = 1/2 
and a = 3/2 , but with the standard definition D"' for the fractional derivative, have 
been discussed by Bagley [30]. Recently, Beyer and Kempfie [46] discussed (4.2) for 
— oo <t< +00 to investigate the uniqueness and causality of the solutions. As they 
let t running in all of IR , they used Fourier transforms and characterized the fractional 
derivative -D" by its properties in frequency space, thereby requiring that for non- 
integer a the principal branch of (iuj)" should be taken. Under the global condition 
that the solution is square summable, they showed that the system described by (4.2) 
is causal iff a > . 

Also here we shall apply the method of Laplace transform to solve the frac- 
tional differential equations and get some insight into their fundamental and impulse- 
response solutions. However, in contrast with the previous section, we now find it 
more convenient to apply directly the formula (1.30) for the Laplace transform of frac- 
tional and integer derivatives, than reduce the equations with the prescribed initial 
conditions as equivalent (fractional) integral equations to be treated by the Laplace 
transform. 
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4.1 The composite fractional relaxation equation 

Let us apply the Laplace transform to the fractional relaxation equation (4.1). 
Using the rule (1.30) we are led to the transformed algebraic equation 

u{s) = co r^ + A^, 0<a<l, 4.3 

wi{s) wi{s) 

where 

wi{s) := s + a s'^ + 1 , (4.4) 

and a > . Putting 

_ 1 + 8°^^^ _ 1 

uo{t) ^ uo{s) := — — , us{t) ^ us{s) := — — , (4.5) 

wi{s) wi{s) 

and recognizing that 

uo(0+) = lim suo{s) = 1 , us{s) = - [suo{s) - 1] , (4.6) 

s— >oo 

we can conclude that 

u{t) = coUo{t)-\- [ q{t-T)us{T)dT, us{t) = - Uo{t) . (4.7) 
Jo 

We thus recognize that uo{t) and us{t) are the fundamental solution and impulse- 
response solution for the equation (4.1), respectively. 

Let us first consider the problem to get 'Uo(t) as the inverse Laplace transform 
of uo{s) . We easily see that the function wi{s) has no zero in the main sheet of the 
Riemann surface including its boundaries on the cut (simply show that Im {wi{s)} 
does not vanish if s is not a real positive number), so that the inversion of the Laplace 
transform uo{s) can be carried out by deforming the original Bromwich path into the 
Hankel path Ha{e) introduced in the previous section, i.e. into the loop constituted 
by a small circle |s| = e with e — > and by the two borders of the cut negative real 
axis. As a consequence we write 



uo{t) = — [ ^ + ' ds. (4.8) 



It is now an exercise in complex analysis to show that the contribution from the 
Hankel path Ha{e) as e — > is provided by 

/■oo 

uo{t)^ / e-^^i7^'J(r;a)dr, (4.9) 
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with 

1 a r" sin (an) 



TT (1 - r)2 + a2 + 2 (1 - r) a cos (ayr) ' 

For o > and < a < 1 the function H^^\{r]a) is positive for all r > since it 
has the sign of the numerator; in fact in (4.10) the denominator is strictly positive 
being equal to |zwi(s)p as s = re^^^ . Hence, the fundamental solution uo{t) has the 
peculiar property to be completely monotone, and H^Q{r;a) is its spectral function. 

Now the determination of us{t) ~ — Wo(t) is straightforward. We see that also the 
impulse-response solution us{t) is completely monotone since it can be represented 

by 

/•oo 

us{t)= e-'"^H^^l,{r;a)dr, (4.11) 

^0 

with spectral function 

(1) , tt{i) /■ \_ 1 ar"sin(a7r) /a to\ 

^a,-i{r; a) = rH^^^ir; = - ^-^ _ + „2 ^2a + 2 (1 - r) ar« cos {an) ' ^^'^^^ 

We recognize that both the solutions uo{t) and us{t) turn out to be strictly 
decreasing from 1 towards as t runs from to 00 . Their behaviour as t — > 0+ and 
t — > 00 can be inspected by means of a proper asymptotic analysis. 

The behaviour of the solutions as t — > 0"'" can be determined from the behaviour 
of their Laplace transforms as Re{s} +00 as well known from the theory of the 
Laplace transform, see e.g. [25]. We obtain as Re{s} +00 , 

uo{s) = s-^ -s-^ + O (s-3+") , us{s) = s-^-a + Q (s-2) , (4.13) 

so that 

-ll — Ct 

uoit) = 1 - t + O {t^-'') , usit) = l-a— r + 0{t), as i^0+. (4.14) 

1 (2 — aj 

The spectral representations (4.9) and (4.11) are suitable to obtain the asymptotic 
behaviour of uo{t) and us{t) as t — > +00, by using the Watson lemma. In fact, 
expanding the spectral functions for small r and taking the dominant term in the 
corresponding asymptotic series, we obtain 

'^o{t)r^a— r, us{t) a— -, as t ^ 00 . (4.15) 

L [1 — a) i {—a) 

We note that the limiting case a = 1 can be easily treated extending the validity 

of eqs (4.3-7) to a = 1 , as it is legitimate. In this case we obtain 

«o(t) = e-V(l + «), ^^(^) = ^e-t/(l + a)^ ^ ^ (4_;Lg) 

1 + a 

Of course, in the case a = we recover the standard solutions uo{t) — us{t) = e~* . 
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We conclude this sub-section with some considerations on the solutions when the 
order a is just a rational number. If we take a = p/q , where p,q & ISl are assumed 
(for convenience) to be relatively prime, a factorization in (4.4) is possible by using 
the procedure indicated by Miller and Ross [10]. In these cases the solutions can be 
expressed in terms of a linear combination of q Mittag-Leffier functions of fractional 
order 1/q, which, on their turn can be expressed in terms of incomplete gamma 
functions, see (A. 14) of the Appendix. 

Here we shall illustrate the factorization in the simplest case a — 1/2 and provide 
the solutions uo{t) and us{t) in terms of the functions ea{t;X) (with a = 1/2), 
introduced in the previous section. In this case, in view of the application to the 
Basset problem, see [24], the equation (4.1) deserves a particular attention. For 
q; = 1/2 we can write 

wi{s) = s + as^/^ + l = {s^/^-X+){s^/^-X-), A± = -a/2±(aV4-l)V2. (4.17) 

Here X± denote the two roots (real or conjugate complex) of the second degree 
polynomial with positive coefficients z'^ + az + 1 , which, in particular, satisfy the 
following binary relations 

A+-A_ = l, A++A_ = -a, A+ -A_ = 2(aV4- 1)^/^ = (a^ -4)1/2. ^^^^g^. 

We recognize that we must treat separately the following two cases 

i) < a < 2 , or a > 2 , and ii) a = 2 , 

which correspond to two distinct roots (A-|- ^ A_), or two coincident roots (A-|- = 
A_ = —1), respectively. For this purpose, using the notation introduced in [24], we 
write 

i) ^ + 



and 



1 + I ^ .V2(,i/2_A+)-.i/2(,i/2_A_)' 

''^ (sV2 + l)2 + 51/2(51/2 + 1)2 ' 

A+ A_ 

^) _1 /9 / _1 /9 \ \ "I" 



where 



^ ^ s + asi/2 + 1 I . 1 ^ 

(5I/2 + 1)2 ' 

A. = ±j^_. (4.21) 



Using (4.18) we note that 

A+ + A_ = l, ^+A_+^_A+ = 0, A+A+ A_ = -a. (4.22) 
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Recalling the Laplace transform pairs (A. 34), (A. 36) and (A. 37) in Appendix, we 
obtain 

i) A_ (A+ Vi) + A+ Ei/2 (A_ Vi) , 



and 



Uo{t) = M{t) := 



U5{t) = N{t) := 



a) (l-2t)Ei/2(-v^) + 2v/tA, 

i) A+ E^/2{X+ Vi) + A_ E^/2{\- Vi) , 
a) (1 + 2t) Ei/2{-Vi) - 2 V^A- 



(4.23) 



(4.24) 



We thus recognize in (4.23-24) the presence of the functions ei/2(t; — A±) = 
^i/2(A± Vi) and ei/2(t) = ei/2(t; 1) = Ei/2i-Vi) ■ 

In particular, the solution of the Basset problem can be easily obtained from (4.7) 
with q{t) = Qo by using (4.23-24) and noting that N{t) dr = 1 - M{t) . Denoting 
this solution by UB{t) we get 



usit) = qo- (go - Co) M{t) . 



(4.25) 



When a = 0, i.e. in the absence of term containing the fractional derivative (due to 
the Basset force), we recover the classical Stokes solution, that we denote by us{t) , 

us{t) = go - {qo - co)e~^ . 

In the particular case qo = cq , we get the steady-state solution UB{t) = us{t) = qo ■ 
For vanishing initial condition cq = , we have the creep-like solutions 

-t] 



UB{t) = qo[l-M{t)] , us{t) = qo 



1 -e 



that we compare in the normalized plots of Fig. 5 of [24] . In this case it is instructive 
to compare the behaviours of the two solutions as t — > 0+ and t — > oo . Recalling the 
general asymptotic expressions of uo{t) = M{t) in (4.14) and (4.15) with a = 1/2, 
we recognize that 



and 



UB{t) = go [t + O {t^^^)] , us{t) = qo [t + O (t^)] , as t ^ 0+ , 



UBit) ^ qo l — a/VTri , us{t) qo — EST] , as t — > oo , 



where EST denotes exponentially small terms. In particular we note that the nor- 
malized plot of UB{t)/qo remains under that of us{t)/qo as t runs from to oo . 

The reader is invited to convince himself of the following fact. In the general case 
< a < 1 the solution u{t) has the particular property of being equal to 1 for all 
t > if g(t) has this property and u{0~^) = 1 , whereas q{t) = 1 for all t > and 
^(0"*") = implies that u{t) is a creep function tending to 1 as t — > oo . 
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4.2 The composite fractional oscillation equation 

Let us now apply the Laplace transform to the fractional oscillation equation 
(4.2). Using the rule (1.30) we are led to the transformed algebraic equations 

a v{s) = co T^ + ci^T + ^Vt' 0<«<1> 4.26a 

W2{S) W2{S) W2{S) 

or 

b v{s)^c^ T^ + ci T^ + ^Vv' l<a<2, 4.266 

W2\S) W2{S) W2{S) 

where 

W2{s) := + a s"^ + 1 , (4.27) 

and a > . Putting 

^o(g):= ^^"^. , 0<a<2, (4.28) 

W2{S) 

we recognize that 

^;o(0+) = lim s^Jo(s) = 1 , = - [svois) - 1] ^ -^o(*) , (4-29) 

s— >oo '?i'2(*j 

and 

= - / voir) dr. 4.30 

«'2(S) S Jo 

Thus we can conclude that 

(a) v{t)^coVo{t)-ciVQ{t)- [ q{t-T)vQ{T)dT, 0<a<l, (4.31a) 

Jo 

or 

(b) v(t) = coVo(t) + ci [ vo{t) dr - [ q{t - t) Vq{t) dr , l<a<2. (4.316) 

In both of the above equations the term ^VQ{t) represents the impulse-response 
solution vs{t) for the composite fractional oscillation equation (4.2), namely the 
particular solution of the inhomogeneous equation with cq = Ci = and with q{t) ~ 
6{t) . For the fundamental solutions of (4.2) we recognize from eqs (4.31) that we have 
two distinct couples of solutions according to the case (a) and (b) which read 

(a) {vo{t),Vie.{t) = -v'o{t)}, (b) {vo{t) , Vi^,{t) ^ [ Vo{T)dT}. (4.32) 

Jo 
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We first consider the particular case a = 1 for which the fundamental and impulse 
response solutions are known in terms of elementary functions. This limiting case 
can also be treated by extending the validity of eqs (4.26a) and (4.31a) to o; = 1 , as 
it is legitimate. From 



vo{s) = 



s + a 



s + a/2 



a/2 



s2 + as+l (s + a/2)2 + (1 - aV4) (s + a/2)2 + (1 - a74) 
we obtain the basic fundamental solution 



(4.33) 



vo{t) 



-at/2 



a 



cos{ujt) H sm{uot) if < a < 2 , 

2uj 



(l-t) if a = 2, 



(4.34) 



,-at/2 



a 



cosh(xt) + — sinh(xt) 



if a > 2 , 



where 



(jj 



(4.35) 



By a differentiation of (4.34) we easily obtain the second fundamental solution via{t) 
and the impulse-response solution vs{t) since vig_{t) = vs{t) = — i'o(^) • We point out 
that all the solutions exhibit an exponential decay as t — > oo . 

Let us now consider the problem to get vo{t) as the inverse Laplace transform of 
vo{s) , as given by (4.26-27), 



vo{t) 



2771 J Br 



sts + as 



a-1 



W2is) 



ds , 



(4.36) 



where Br denotes the usual Bromwich path. Using a result by Beyer and Kempfle 
[46] we know that the function W2 {s) (for a > and 0<q:<2,q:^1) has exactly 
two simple, conjugate complex zeros on the principal branch in the open left half- 
plane, cut along the negative real axis, say s+ = pe and s_ = pe~'^'^ with p > 
and 7r/2 < 7 < TT. This enables us to repeat the considerations carried out for the 
simple fractional oscillation equation to decompose the basic fundamental solution 
vo{t) into two parts according to vo{t) = fa{t; a) + Qait; a) . In fact, the evaluation 
of the Bromwich integral (4.36) can be achieved by adding the contribution fa{t] a) 
from the Hankel path Ha{e) as e — > , to the residual contribution ga{t; a) from the 
two poles s± . 

As an exercise in complex analysis we obtain 



/■oo 

fa{t;a)= / e-^*iy('J(r;a)dr 

^0 



(4.37) 
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with spectral function 

TT I _ I 

(4.38) 



HaQ{r;a) = - -Im i -— ^ 

TT t Ms) ,=^e-J 

1 a r""-*^ sin (ctTr) 



TT (r2 + 1)2 + a2 r^a + 2 (r^ + 1) a r« cos (cktt) ' 

Since in (4.38) the denominator is strictly positive being equal to |w2(-s)P as s = 
re , the spectral function H^Q{r;a) turns out to be positive for all r > for 
< CK < 1 and negative for all r > for 1 < o; < 2 . Hence, in case (a) the function 
fait) , in case (b) the function —fait) is completely monotone; in both cases fait) 
tends to zero as t — > cxd , from above in case (a), from below in case (b), according to 
the asymptotic behaviour 

fait; a) ~ a — — r , as t ^ oo , 0<q:<1, 1<q;<2, (4.39) 

r(l — a) 

as derived by applying the Watson lemma in (4.37) and considering (4.38). 
The other part, gait; a) , is obtained as 

gait; a)=e^+'^Res 



s + as"^ 
W2is) 



+ conjugate complex 



(4.40) 



2 s+ + a a ^ 



Thus this term exhibits an oscillatory character with exponentially decreasing am- 
plitude like exp (— | cos 7I) . 

Then we recognize that the basic fundamental solution fo(t) exhibits a finite 
number of zeros and that, for sufficiently large t , it turns out to be permanently 
positive if < a < 1 and permanently negative if 1 < a < 2 with an algebraic decay 
provided by (4.39). 

For the second fundamental solutions via(i) , ■i'ib(*) and for the impulse- response 
solution vsit) , the corresponding analysis is straightforward in view of their connec- 
tion with Voit), pointed out in (4.31-32). The algebraic decay of all the solutions as 
t— >oo,forO<Q;<l and 1 < ck < 2 , is henceforth resumed in the relations 

-f-—ce ^ — -j-l — a 

^o(*)~a— -, vi^it) = vsit) ^ -a— -, vii,it)r^a— -. (4.41) 

L [1 — a) i {—a) i (2 — q;) 

In conclusion, except in the particular case ct = 1 , all the present solutions of 
the composite fractional oscillation equation exhibit similar characteristics with the 
corresponding solutions of the simple fractional oscillation equation, namely a finite 
number of damped oscillations followed by a monotonic algebraic decay as t — > 00 . 
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5. CONCLUSIONS 

Starting from the classical Riemann-Liouville definitions of the fractional integra- 
tion operator and its left-inverse, the fractional diff'erentiation operator, and using the 
powerful tool of the Laplace transform method, we have described the basic analytical 
theory of fractional integral and differential equations. For a numerical treatment we 
refer to Gorenflo [23] and to the references there quoted. 

For the fractional integral equations we have considered the basic examples pro- 
vided by the linear Abel equations of first and second kind. For both the kinds we 
have given the solution in different forms and discussed an interesting application to 
inverse heat conduction problems. 

Then we have analyzed in detail the scale of fractional ordinary differential equa- 
tions (FODE), see (3.5), u{t) + u{t) = q{t), t>0, < a < 2 , with a modified 
fractional differentiation , the Caputo fractional derivative, that takes account of 
given initial values tt(0"'") if < a < 1 , the case of fractional relaxation, tt(0"'") and 
u'{Q'^) if 1 < a < 2 , the case of fractional oscillation. 

We have investigated in depth the properties of the fundamental and the impulse- 
response solutions. All these solutions can be explicitly written down in terms of 
Mittag-Leffter functions. They tend to zero like powers t~^ (with appropriate choices 
of P) , monotonically if < a < 1 , but exhibiting finitely many oscillations around 
zero if 1 < q: < 2 (the more of these the nearer a is to the limiting value 2). 
If 1 < a < 2 these equations are able to model processes intermediate between 
exponential decay {a = 1) and pure sinusoidal oscillation (a = 2). We have found 
these qualitative properties essentially by bending the Bromwich integration path of 
the Laplace inversion formula into the Hankel path, thus for each of these functions 
obtaining an integral representation as the Laplace transform of a function that 
nowhere changes its sign, augmented if 1 < a < 2 by an oscillatory contribution 
resulting from a pair of conjugate complex poles lying in the left half-plane. 

By quite analogous methods we have studied the composite equations, see (4.1-2), 
{D + aD^ + 1) u{t) = q{t) , < a < 1 , and (I^^ ^ aD^ + 1) v{t) = q{t) , < a < 2 , 
with a > , which model processes of relaxation and of oscillation, respectively. We 
have obtained similar properties of the fundamental and impulse-response solutions 
with respect to monotonicity and oscillatory behaviour. 

Let us stress the fact that our adoption of the Caputo fractional derivative D'^ 
with m — 1 < a < m , meN, and the consequent prescription of the initial values 
in analogy with the ordinary differential equations of integer order m , stands in 
contrast to the majority of the treatments of fractional differential equations, where 
the standard fractional derivative is used, see e.g. [5], [10]. As pointed in §1.3 the 
adoption of D" requires the prescription of certain fractional integrals as t — > 0"*" . 
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In our opinion, the different prescription of the the initial data points out the 
major difference between the two definitions for the fractional derivative. The analogy 
with the cases of integer order would induce one to adopt the Caputo derivative in 
the treatment of differential equations of fractional order for physical applications. In 
fact, in physical problems, the initial conditions are usually expressed in terms of a 
given number of bounded values assumed by the field variable and its derivatives of 
integer order, no matter if the governing evolution equation may be a generic integro- 
differential equation and therefore, in particular, a fractional differential equation. 

The liveliness of the field of fractional integral and differential equations, both in 
applications and in pure theory, is underlined by several papers and some books that 
have appeared recently. We would like to conclude the present lectures with brief 
hints to recent investigations, which have not been quoted explicitly up to now since 
not strictly related to our results, but have attracted our attention. Naturally, this 
listing is far from exhaustive. The interested reader can find more on problems and 
aspects in several papers recently published or in press in some conference volumes 
and specialized magazines. 

We first like to quote the most recent book by Rubin [14], who starting from 
one-dimensional fractional calculus develops the theory of multidimensional weakly 
singular integral equations of first kind, making heavy use of the Marchaud approach. 
We thus recognize that all existing books on fractional calculus vary widely from each 
other in their character concerning problems treated and methods applied. We quote 
also the Ph.D. thesis by Michalski [47], who treats linear and nonlinear problems of 
fractional calculus (in one and in several dimensions) in a very elegant way. 

The importance of using fractional methods in physics for describing slow decay 
processes and processes intermediate between relaxation and oscillation was stressed 
by Nigmatullin [48] in 1984. Nonnenmacher and associates published a series of 
papers (of which we quote [49-50]) discussing various physical aspects of fractional 
relaxation. Fractional relaxation is overall a peculiarity of a class of viscoelastic bodies 
which are extensively treated by Mainardi [24] , to which we refer for details and addi- 
tional bibliography. The fractional calculus finds important applications in different 
areas of applied science including electrochemistry, see e.g. [51-54], electromagnetism, 
see e.g. [55-56], radiation physics, see e.g. [57-59], and control theory, see e.g. [60-64]. 

Yet another field of applications of fractional calculus is that of fractional partial 
differential equations (FPDE), including certain equations of fractional diffusion, 
introduced to explain the phenomena of anomalous diffusion in complex or fractal 
systems. We refer again to Mainardi [24] for a mathematical treatment of a rele- 
vant FPDE, referred to as the time fractional diffusion-wave equation, with some 
applications and related references. 
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APPENDIX: THE MITTAG-LEFFLER TYPE FUNCTIONS 

In this Appendix we shall consider the Mittag-Leffler function and some of the 
related functions which are relevant for their connection with fractional calculus. It 
is our purpose to provide a review of the main properties of these functions including 
their Laplace transforms. 

A.l The Mittag-Leffler functions Ea{z) , Ef^^p^z) 

The Mittag-Leffler function Ea{z) with a > is so named from the great Swedish 
mathematician who introduced it at the beginning of this century in a sequence of 
five notes, see [65-69] . The function is defined by the following series representation, 
valid in the whole complex plane. 



n=0 ^ 

It turns out that Ea{z) is an entire function of order p = 1/a and type 1. This 
property is still valid but with p = 1 /Re{Q!} , if o: e C with positive real part, as 
formerly noted by Mittag-Leffler himself in [68] . 

In the limit for a ^ the analyticity in the whole complex plane is lost since 

oo ^ 
n=0 

The Mittag-Leffler function provides a simple generalization of the exponential 
function because of the substitution of n\ = T{n + 1) with (an)! = T{an + 1) . 
Particular cases of (A.l), from which elementary functions are recovered, are 

E2 {+z^) = cosh z , E2 {-z^) = cos 2; , z eV, {A.3) 

and 

Ei/2(±^'/2) = e^ [1 + erf (i^'/^)] = e^ erfc (tz^/^) , zeC, (AA) 
where erf (erfc) denotes the (complementary) error function defined as 

erf (z) := ^ / e""^ du , erfc {z) := 1 - erf (2) , zeC. 

V^r Jo 

In (A. 4) for z^^'^ we mean the principal value of the square root of z in the complex 
plane cut along the the negative real axis. With this choice ±z^/'^ turns out to be 
positive/negative for z e R"*". 

Since the identities in (A.3) are trivial, we present the proof only for (A. 4). Avoid- 
ing the inessential polidromy with the substitution ±2;^/^ — > 2; , we write 
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Whereas the even part is easily recognized to be u{z) = exp{z^) , only after some 
manipulations the odd part can be proved to be v{z) = exp(2;^) erf (2;) . To this end 
we need to recall the following series representation for the error function, see e.g. 
the handbook of the Bateman Project [70] or that by Abramowitz and Stegun [71] , 



erf(^) = ^e-^ V — ^— — + 1 
VtF ^„ 2m+l !! 

and note that (2m + 1)!! := 1 • 3 • 5 . . . • (2m + 1) = 2"*+^ r(m + 3/2)70?. An 
alternative proof is obtained by recognizing, after a term-wise differentiation of the 
series representation in (A. 5), that v{z) satisfies the differential equation in C , 

^ ^zv{z) , ^;(0) = 0, 



whose solution can immediately be checked to be 

2 2 /"^ _ 2 2 
v{z) = —=Q^ \ e du = e^ erf (2;). 
Jo 

A straightforward generalization of the Mittag-Leffler function, originally due to 
Agarwal in 1953 based on a note by Humbert, see [72-74], is obtained by replacing the 
additive constant 1 in the argument of the Gamma function in (A.l) by an arbitrary 
complex parameter /? . Later, when we shall deal with Laplace transform pairs, the 
parameter /3 is required to be positive as a . For the new function we agree to use 
the following notation 

00 „ 

EaA^) ■■= E r(Z ^ m ' «>0'/?eC, zeC. (A6) 

n=0 ^ 



Particular simple cases are 



e 



z 



1 ^ , , sinh(2i/2) 



Ei,2{z) = ^—, E2,2{z) = ^7^. (A7) 

We note that E^^z) is still an entire function of order p = 1/a and type 1 . 

In these lectures we have preferred to use only the original Mittag-Leffler function 
(A.l) since our problems depend on only a single parameter a, the order of fractional 
integration of differentiation. However, for completeness, we list hereafter the general 
functional relations for the generalized Mittag-Leffler function (A. 6), which involve 
both the two parameters a , /? , see [18] and [70], 

EaA^) = + ^ ^a,/3+a(^) , (^-8) 

EaA^) = f3Ea,f3+l{z) + «^ ^ Ea,p+liz) , (A.9) 

(J^J [zf'-'E^An] = z^-^-^E^,p_^{z'^) , p e N. (^.10) 
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A 2. The Mittag-Leffler functions of rational order 

Let us now consider the Mittag-Leffler functions of rational order a = p/q with 
p, g e ]N relatively prime. The relevant functional relations, that we quote from 
[18], [70], turn out to be 



dzP 



Ep/.izP/'^ 



q-l 



^-kp/q 



r(l - kp/q) 



g = 2,3. 



p-i 



and 



h=0 

q-l 



^ ^{l-k/q,z) 
+ T{l-k/q) 



g = 2,3,. 



(All) 
(A12) 

(A13) 
{A.U) 



where 7(0, z) := Jq e u°'~^ du denotes the incomplete gamma function. Let us now 
sketch the proof for the above functional relations. 

One easily recognizes that the relations (A. 11) and (A. 12) are immediate conse- 
quences of the definition (A.l). 

In order to prove the relation (A. 13) we need to recall the identity 



p-i 

h=0 



^i2nhk/p _ 



p if A; = (mod p) , 
if ^ (mod p) . 

In fact, using this identity and the definition (A.l), we have 

p-i 

J2Ea{ze'^^^/n=pEap{zn, 

Substituting in the above relation a/p instead of a and z^^^ instead of 2; , we obtain 

1 ^"^ 



(A15) 



(A16) 



^=0 



Setting above a =p/q, we finally obtain (A. 13). 

To prove the relation (A. 14) we consider (A. 12) for p = 1 . Multiplying both sides 
by , we obtain 

.-k/q 



d_ 

dz 



El/q[z 



— e 



^ rfi 



(^.18) 



.=1 r(i - k/,) ■ 

Then, upon integration of this and recalling the definition of the incomplete gamma 
function, we arrive at (A. 14). 



266 Fractional Calculus: Integral and Differential Equations of Fractional Order 



The relation (A. 14) shows how the Mittag-Leffler functions of rational order can 
be expressed in terms of exponentials and incomplete gamma functions. In particular, 
taking in (A. 14) q = 2 , we now can verify again the relation (A. 4). In fact, from 
(A. 14) we obtain 

1 + ^7(1/2,2;) , (A19) 

which is equivalent to (A. 4) if we use the relation erf (z) = ^{1/2, z^)/-\/n , see e.g. 
[70-71]. 



Ey.iz'/') = 



A3. The Mittag-Leffier integral representation and asymptotic expansions 

Many of the most important properties of E^iz) follow from Mittag-Leffler's 
integral representation 

Ea{z)^^ [ ^y^^dC, a>0, ^eC, (A20) 

^TTZ J Ha C - ^ 

where the path of integration Ha (the Hankel path) is a loop which starts and ends at 
—00 and encircles the circular disk |C| < \z\^/°' in the positive sense: — tt < arg^ < tt 
on Ha . To prove (A. 20), expand the integrand in powers of integrate term-by-term, 
and use Hankel's integral for the reciprocal of the Gamma function. 

The integrand in (A. 20) has a branch-point at C = 0. The complex ^-plane is 
cut along the negative real axis, and in the cut plane the integrand is single- valued: 
the principal branch of is taken in the cut plane. The integrand has poles at the 
points — -2^^" Q^TTtm/a ^ ^ integer, but only those of the poles lie in the cut plane 
for which — cktt < arg 2 + 27rm < an . Thus, the number of the poles inside Ha is 
either [a] or [a -|- 1] , according to the value of arg z. 

The integral representation of the generalized Mittag-Leffler function turns out 
to be 

1 r ^"-/SgC 



E^,f3{z) = — / ^ dC, «,/3>0, zee. (A21) 

^^'i J Ha S ~ Z 

The most interesting properties of the Mittag-Leffler function are associated with 
its asymptotic developments as ^ ^ cxd in various sectors of the complex plane. These 
properties can be summarized as follows. 
For the case < ct < 2 we have 

£;«(2;)~-exp(2;i/«)-^— ^— — , |2;|^oo, | arg 2;] < a7r/2 , (A22) 

00 _ t. 

Ea{z)^-'^Y(l IT' \z\^oo, a7r/2 < arg 2; < 27r - Q;7r/2 . {A.23) 

k=i ^ " 
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For the case a > 2 we have 

Eaiz) ~ - ^exp ^V«e2--/a] _ ^ , \z\ ^ oc , iA.24) 

m fc=l ^ ' 

where m takes all integer values such that —a'K/2 < arg 2; + 27rm < a7r/2 , and arg z 
can assume any value between — tt and +7v inclusive. 

Prom the asymptotic properties (A. 22-24) and the definition of the order of an 
entire function, we infer that the Mittag-Leffler function is an entire function of order 
1/a for a > 0; in a certain sense each Ect{z) is the simplest entire function of its 
order, see Phragmen [75]. The Mittag-Leffler function also furnishes examples and 
counter-examples for the growth and other properties of entire functions of finite 
order, see Buhl [76]. 

A4- The Laplace transform pairs related to the Mittag-Leffler functions 

The Mittag-Leffler functions are connected to the Laplace integral through the 
equation 

rOO 1 rOO 

/ e-'^ Ea {u"^ z) du = = / e-" w^"^ E^piu'^z) du, a, P>0. {A.25) 

Jo 1 - ^ Jo 

The integral at the L.H.S. was evaluated by Mittag-Leffler who showed that the region 
of its convergence contains the unit circle and is bounded by the line Kez^^"^ — 1. 

The above integral is fundamental in the evaluation of the Laplace transform of 
-Eq, (—At") and -Ea./J (^-^^") with a , /3 > and A G C . Since these functions turn 
out to play a key role in problems of fractional calculus, we shall introduce a special 
notation for them. 

Putting in (A.25) u — st and u!^ z = —At" with t > and A e C , and using the 
sign for the juxtaposition of a function depending on t with its Laplace transform 
depending on s, we get the following Laplace transform pairs 

e«(t; A) := E^ (-A^) - , Res> |A| V« , (^.26) 

S -r A 

and ^ 

ea,/3(t; A) := tl"-^ Ea,^ (-AT) - , Res> \X\^/^ . {A.27) 

We note that the results (A. 26-27), but with a different notation, were used by 
Humbert and Agarwal [72-74] to obtain a number of functional relations satisfied by 
Ea{z) and Eajj{z) . Of course the results (A. 26-27) can also be obtained formally 
by Laplace transforming term by term the series (A.l) and (A. 6) with z = —At" , 
respectively, and summing the resulting series. 
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We find worthwhile to hst the following relations for the functions Ca,/} easily- 
obtained from (A. 8-9) : 

ea,f3it;X) = - Xea,(3+ait;X) , {A.28) 

and 

— ea,/3+i(t;A) = ea,^(t;A). (A29) 

A remarkable property satisfied by the functions ea{t; X) , Ca^/sit; X) when A is 
positive and 0<q;<1, < a < (3 < 1 , respectively, is to be completely monotone 
for t > . 

We recall that a function f{t) is told to be completely monotone for t > if 
(t) > for all n = 0,1,2,... and for all t > , and that a sufficient 
condition for this is the existence of a nonnegative locally integrable function K{r) , 
r > , referred to as the spectral function, with which f{t) = e K{r) dr . For 
more details see e.g. the book by Berg & Forst [77]. 

Excluding the trivial case a = /5 = 1 for which ei(t; A) = ei,i(t; A) = e"'^* , we 
can prove the existence of the corresponding spectral functions using the complex 
Bromwich formula to invert the Laplace transform in (A. 26-27) and bending the 
Bromwich path into the Hankel path, as we have already shown in the special case 
ea{t) :— Cait] 1) in §3. As an exercise in complex analysis (we kindly invite the reader 
to carry it out) we obtain the integral representations [A. 30-33], 

e«(t;A):=/ e~^^ Ka{r; X) dr , 0<a<l, A>0, (A30) 
Jo 

with spectral function 

Ka(r; A) = - — \ — ( — > , (A31) 

«v > ) TT r2« + 2 A r« cos (avr) + A2 - ' ^ ' 

and 

/•oo 

ea,/3(^;A):= / e-^^is:«,/3(r;A)dr, 0<a</3<l, A>0, (A32) 
Jo 

with spectral function 

,{r; A) = i A sin [(/? - «).] + r" sin (/j.) ^ « - > q . ^^ 33) 

' ' ^ r2« + 2Ar" cos(Q;7r) + A2 - ^ ' 
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Historically, the complete monotonicity of the Mittag-Leffler function in the neg- 
ative real axis, i.e. of Ea{—x) , for a; e H"*", when < a < 1 , was first conjectured 
by Feller using probabilistic methods and rigorously proved by Pollard in 1948 [78]. 

Only recently, Schneider [79] has proved a theorem for the complete monotonicity 
of the generalized Mittag-Leffler function in the negative real axis. He proved that 
Ea,f3{—x) , for X e R"*", is completely monotone iffO<a<l and P > a. Our 
conditions for ea,/3(t, A) to be completely monotone appear more restrictive than 
those by Schneider for Eo,^f3{—x) ; however, we must note that in our case (A. 27) the 
factor t^~^ precedes the generalized Mittag-Leffler function. 

We note that, up to our knowledge, in the handbooks containing tables for the 
Laplace transforms, the Mittag-Leffler function is ignored so that the transform pairs 
(A. 26-27) do not appear if not in the special cases a = 1/2 and /? = 1 , 1/2 , written 
however in terms of the error and complementary error functions, see e.g. [71]. In 
fact, in these cases we can use (A. 4) and (A. 28) and recover from (A. 26-27) the two 
Laplace transform pairs 

f ei/2(t;±A) =e'^^^erfc(±AVt), AeC, {A.M) 

(^.35) 



sV2 (si/2 ± X) 

1 . . 1 



1/2 , X ^ ei/2,1/2 {t; ±A) = ^= T A ei/2 (t; ±A) , AeC. 
We also obtain the related pairs 



sV2 (sl/2 ± A)2 



^ 2y^T2Atei/2(t;±A), AeC, 



(^.36) 



^^i/2^;,)2 - T2A-^- + (l + 2A^t)ei/2(t;±A), AeC, (A37) 
In the pair (A. 36) we have used the properties 



si/2(si/2±A)2 ds\sy^±Xj' ds 

The pair (A. 37) is easily obtained by noting that 



1 1 A 



(sV2 ± A)2 si/2 (si/2 ± A) ^ sV2 (si/2 ± xy 
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A. 5 Additional references for the Mittag-Leffier type functions 

We note that the Mittag-Leffler type functions are unknown to the majority of sci- 
entists, because they are ignored in the common books on special functions. Thanks 
to our suggestion the new 2000 Mathematics Subject Classification has included these 
functions, see the item 33E12: Mittag-Leffier functions and generalizations. 

A description of the most important properties of these functions with relevant 
references can be found in the third volume of the Bateman Project [70], in the chap- 
ter XVIII devoted to miscellaneous functions. The specialized treatises where great 
attention is devoted to the Mittag-Leffier type functions are those by Dzherbashyan 
[18], [22]. For the interested readers we also recommend the classical treatise on com- 
plex functions by Sansone & Gerretsen [80] , where a sufficiently detailed treatment of 
the original Mittag-Leffler function is given. Since the times of Mittag-Leffler several 
scientists have recognized the importance of the Mittag-Leffler type functions, pro- 
viding interesting results and applications, which unfortunately are not much known. 
As pioneering works of mathematical nature in the field of fractional integral and 
differential equations, we like to quote those by Hille & Tamarkin and by Barret. In 
1930 Hillc & Tamarkin [81] have provided the solution of the Abel integral equation 
of the second kind in terms of a Mittag-Leffler function, whereas in 1956 Barret [82] 
has expressed the general solution of the linear fractional differential equation with 
constant coefflcients in terms of Mittag-Leffler functions. 

As former applications in physics we like to quote the contributions by K.S. Cole 
(1933), quoted by H.T. Davis [15, p. 287] in connection with nerve conduction, and 
by F.M. de Oliveira Castro (1939) [83] and B. Gross (1947) [84] in connection with 
dielectrical and mechanical relaxation, respectively. Subsequently, in 1971, Caputo 
& Mainardi [28] have proved that the Mittag-Leffler function is present whenever 
derivatives of fractional order are introduced in the constitutive equations of a linear 
viscoelastic body. Since then, several other authors have pointed out the relevance 
of the Mittag-Leffler function for fractional viscoelastic models, see Mainardi [24]. 

In recent times the attention of mathematicians towards the Mittag-Leffler type 
functions has increased from both the analytical and numerical point of view, overall 
because of their relation with the fractional calculus. In addition to the books and 
papers already quoted in the text, here we would like to draw the reader's attention 
to the most recent papers on the Mittag-Leffler type functions, e.g. Al Saqabi & 
Tuan [85], Kilbas &; Saigo [86], Gorenflo, Luchko & Rogozin [87] and Mainardi &; 
Gorenflo [88] . Since the fractional calculus has actually recalled a wide interest for 
its applications in different areas of physics and engineering, we expect that soon the 
Mittag-Leffler function will exit from its isolated life as Cinderella (using the term 
coined by F.G. Tricomi in the 1950s for the incomplete Gamma function). 
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